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FINITE ELEMENT ANALYSIS OF TRANSONIC FLOWS IN 
CASCADES - IMPORTANCE OF COMPUTATIONAL GRIDS 
IN IMPROVING ACCURACY AND CONVERGENCE 


Akin Ecer* and Hasan U. Akay** 

Purdue University at Indianapolis 

SUMMARY 

The finite element method is applied for the solution of transonic 
potential flows through a cascade of airfoils. Convergence characteristics 
of the solution scheme are discussed. Accuracy of the obtained numerical 
solutions is studied for various flow regions in the transonic flow con- 
figuration. The design of an efficient finite element computational grid 
is discussed in improving accuracy and convergence. 

This study shows that the solution of nonlinear equations corresponds 
to the numerical integration of a pseudo- time system in the solution of 
steady, inviscid potential equations. The choice of the artificial vis- 
cosity parameter, which is necessary for obtaining convergence in the super- 
sonic pocket, and the relaxation parameter, which controls both the conver- 
gence characteristics and the accuracy of the obtained solution, are inves- 
tigated. Proper choice of these parameters as a function of the computational 
grid is presented. The accuracy of the obtained solutions around the leading 
and trailing edges of the airfoil, in the subsonic and supersonic flow regions, 
and around the shock are discussed, and appropriate grid refinement techniques 
are suggested. An efficient numerical integration scheme for the solution of 
transonic flow equations is also presented. 

The numerical results include several basic test cases to demonstrate 
the degree of grid refinement required for accuracy and convergence. Also, 
included in the report are results for complicated flow configurations with 
choked flow conditions and cascades with high staggers and high solidities. 


* Professor of Mechanical Engineering 

** Research Associate 



INTRODUCTION 


Analysis of transonic flows through a cascade of airfoils requires 
accurate modeling of the flow boundaries as well as the details of the 
subsonic and supersonic flow regions and the shock. The following basic 
properties must be considered in the analysis of this problem: 

. Equations change their characteristics in 
subsonic and supersonic flow regions which 
have irregular shapes and may be separated 
by a discontinuity, 

. Leading and trailing edges of the airfoils 
are sources of near singularities in the flow, 

. Cascades usually have irregular geometries . 

The efficiency of the computational techniques in treating such irregular 
geometric configurations is directly related to the applicability of 
irregular computational grids. When such grids are designed properly, 
sufficiently accurate approximations of such irregular flow configurations 
can be obtained. 

Numerical solutions of steady, inviscid transonic flow problems have 
been treated using the finite difference method by several researchers 
during the last decade (refs. 24,25,28,29). More recent studies have been 
in the direction of improving the efficiency and accuracy of these solutions 
(refs. 8,10,16,19,20,21,23,26,31,32). During these investigations, the 
importance of a computational grid in the calculation of transonic flows 
was observed. The modification of computational grids has been explored 
as one of the main tools for improving efficiency and accuracy (refs. 21, 
26,31). 

Application of the finite element method for the analysis of inviscid 
transonic flow problems has also been considered by several investigators 
during the last five years, e.g., (refs. 1-4,6,7,9,12,15,17). The impor- 
tance of the computational grid in obtaining an accurate solution with 
finite element method was reported by the present authors in reference 15. 
One of the main advantages of the finite element method has been its ability 
to accommodate irregular grids. A less utilized advantage is the capability 
of adapting the computational grid to the solution either automatically or 
by interactive means. This property becomes quite important in solving 
transonic flows in cascades and in choosing the most suitable computational 
grid for a particular flow problem. 

The objectives of the present study are to develop efficient computa- 
tional techniques for the solution of transonic flows in cascades, and to 
investigate the important features of computational grids for this problem. 
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The existing literature on the development of a computational grid for 
improving accuracy and convergence of numerical solutions has dealt 
basically with elliptic systems. Accuracy considerations are generally 
based on refining regions with steep gradients. For improving conver- 
gence characteristics of nonlinear elliptic equations, multigrid methods 
have been developed recently to smooth the residiials with high frequency 
content in relaxation procedures (ref. 5) . 

Due to the hyperbolic nature of the equations in the supersonic 
pocket, the development of an efficient computational grid for analyzing 
transonic flow through a cascade of airfoils cannot be based only on the 
above considerations. In the computations it is observed that the con- 
vergence problem is far more critical in the supersonic pocket than it 
is in the surrounding subsonic flow region. Also, convergence and accuracy 
are generally most critical around the shock. 

In this report, the basic considerations in the design of a finite 
element grid for obtaining accurate and efficient solutions to the tran- 
sonic flow problem for cascades are presented. Effects of the grid on 
subsonic and supersonic flow regions and the shock are investigated sepa- 
rately. For an accurate and efficient solution, a computational procedure 
should be based on consideration of each of these problems and their 
interaction . 

In the following chapter, the finite element analysis of the transonic 
flow problem for a cascade of airfoils is summarized. In particular, the 
treatment of periodic boundary conditions and the Kutta condition for sharp 
trailing edges are discussed. A combined shock fitting and shock capturing 
scheme is presented. It is also shown that the solution of the set of non- 
linear equations for steady transonic flows requires a numerical integration 
in time of a pseudo-time dependent system of equations. A constant coef- 
ficient integration scheme is developed to improve the computational 
efficiency of this numerical integration scheme. 

In the next chapter, convergence characteristics to a steady state are 
discussed. For the subsonic region, it is shown that the scheme basically 
requires the solution of an elliptic equation with a relaxation procedure. 
For the supersonic flow region, the convergence is both a function of the 
relaxation parameter w and the artificial viscosity parameter which is 

introduced by upwinding the element density. It is emphasized that while 
convergence improves with increasing values of a^, the accuracy of the 

solutions is very sensitive to the amount of artificial viscosity. The 
relationship between the computational grid and the distribution of arti- 
ficial viscosity is found to be very important and is discussed in detail. 
This convergence study is also applicable to finite difference and finite 
volume methods based on the upwinding of local density. 

In the next chapter, basic concepts in the development of computational 
grids for the solution of transonic flow problems are related to the ac- 
curacy of the numerical solutions. For obtaining accurate solutions in 
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subsonic flow regions, grid refinement of regions with high pressure 
gradients is important. In the supersonic pockets, where pressure gradients 
are generally small, the accuracy problem is mostly related to the appli- 
cation of an artificial viscosity. Simple yet efficient finite element 
grids can be employed for the solution of this problem if the artificial 
viscosity is correctly applied. It is shown that accuracy around the shock 
is also very sensitive to the application of the artificial viscosity. 
Finally, the sensitivity of the flow at the trailing edge is discussed 
for cases where a shock is close to the trailing edge. 

The last chapter includes several case studies illustrating the 
generality of the developed computational scheme and the effect of com- 
putational grids in terms of accuracy and convergence. The test cases 
include a) transonic flow around an isolated cylinder, b) a cascade of 
NACA 0012 airfoils, with a stagger angle of 45° and gap to chord ratio 
of 3.6, c) a cascade of NACA 0012 airfoils, with a stagger angle of 45° 
and gap to chord ratio of 1.0, (choked flow), d) a cascade of 8% 
circular-arc airfoil, with a stagger angle of 45° and gap to chord ratio 
of 1.0, and e) a cascade of Jose Sanz stator blades. 

The work presented in this report summarized the applicability of 
the finite element method for the solution of transonic flow in cascades. 

The ability to interactively modify the finite element grid in order to 
improve the accuracy in specific flow regions proves to be one of the main 
advantages of the method. It is demonstrated in this report that a 
designer can analyze transonic flows in cascades: 

. by starting with a rough finite element grid and 
refining the critical regions if necessary, 

. by controlling the artificial viscosity distri- 
bution over the grid and through the iterations, 

. by using a combined shock fitting and shock capturing 
scheme to check the accuracy of the shock region, 

. by checking the accuracy of the solution from the 
comparison of results obtained using different grids. 

By using the above procedure, one can always ensure the correctness of the 
obtained results. It should also be noted that, as the problems get more 
complex, the need for such a programmed approach becomes more important. 
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FINITE ELEMENT FORMULATION OF TRANSONIC FLOWS 


Governing Equations 

We consider the solution of the full potential equation for steady, 
inviscid and irrotational flows expressed in the conservative form 

(p(|) ) + (p 4) ) = 0 in - (1) 

where (j) is the velocity potential, Q is the solution domain, and p is the 
mass density of the fluid. The general boundary conditions of equation (1) 
are 


(p = (p^ on Si (2) 

and 

pcj) = f on S 2 , (3) 

where and f are some specified quantities on the boundary S = Si + S 2 , 
and n is the outward normal on S 2 . 

By combining the isentropic equation of state with the Bernoulli equa- 
tion, one obtains the relation 

1 

p = const (K^ - q^) ^ ^ (4) 

for the mass density, where y is the ratio of specific heats. 


and 




2a^ 

Y-1 


+ 




(5) 

( 6 ) 


In the above, 
out the flow. 


K is the maximum attainable speed which is constant through- 
a is the local speed of sound, and q is the local flow speed. 
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For the cascade geometries shown 

sound a. is normalized then equation 
in 


in figure 1, 
(5) becomes 


if the inlet speed of 


= -At + M . ^ 

y-1 in 


where M, = q. /a. is the inlet Mach number, 
in in in 


( 7 ) 


Cascade Geometry 

Either of the solution domains shown in figure 1 can be used to model 
two-dimensional cascades. In each case, a set of periodic boundaries to- 
gether with the blade surfaces, the inlet and exit boundaries, the stagger 
angle and the pitch h define the cascade geometry. Since the domain in 

figure lb is multiply connected, a splitting boundary such as the slit line 
E-F shown in the figure is required in order to allow the velocity potential 
to be discontinuous across E-F. This assures that (}) will be single valued. 

The configuration shown in figure lb is found to be more advantageous 
for highly staggered cascades. Firstly, it lends itself easily to the 
generation of less distorted grids, especially at the rounded leading and 
training edges. Secondly, since the periodic boundaries are furthest from 
the blade surfaces, the effects of approximations introduced in satisfying 
the periodicity conditions in this case are minimal. 


The boundary conditions for configuration of figure la have been dis- 
cussed in reference 15. Only the case b will be discussed here. Assuming 

that the farfield inlet and exit flow conditions are uniform and M. , 3. , 

in in 

3g^ are known, the boundary conditions become 


i) 

f = 

f . = 
in 

-p , q , CosB . on A-B 

in in in 

(8) 

ii) 

f = 

f 

ex 

- f , on C-D 

in 

(9) 

iii) 

f = 

0 on 

the airfoil surface 

(10) 

iv) 

- 

■ ({) = 

X on E-F 

(11) 


where X is 
ditions as 

X 


the 

= q 


circulation which may 

. (h Cos3. + h Sin3 . ) 
in X in y in 


be 


determined from the farfield con- 

q (h Cos3 + h Sin3 ) (12) 

ex x ex y ex 
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where 


h = hSina , h = hCosa , 
X s . y s 


( 13 ) 


and is computed iteratively from the continuity equation 

f. + f =0. 
in ex 

The periodicity condition is: 


v) (J)(x+h^,y+h^) 


= (j) (x,y) + q, (h Cos3 . + h Sin3. ) on B-D (15) 

1 in X in y in 

where (p (x,y) is the unknown distribution of (p on A-C. Finally, 

1 

vi) (|) = 0 on A . (16) 

In the above, equation (9) is due to conservation of mass, (15) is due to 
perodicity, and (16) is needed to avoid the arbitrariness of the potential 
function (p. 

In some cases the exit flow angle 3 is not known a priori; therefore, 

ex 

the circulation A is also an unknown. To determine A in this case, the 
Kutta condition 


q^l + = q"l . (17) 

E E 

must be satisfied at the sharp trailing edge of the airfoil. A convenient 
way of satisfying equation (17) and determining A in an iterative fashion 
is discussed on page 18. 


Finite Element Formulation 


For a finite element formulation of equation (1) - (3) we use the 
Bateman’s variational functional (ref. 13) 

TT = f^p d^2 + f(p dS 
b2 


(18) 


where for isentropic and perfect gases 

- 1=1 


P = 2 ^ const p 


(19) 


is the pressure, 
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It can be verified that Ritz type finite elements obtained through the 
consistent use of Bateman’s functional are equivalent to applying a Galerkin 
formulation to equation (1) (ref. 1). 

To introduce the finite element approximations, we express equation (18) 
for an element e and sum the contributions of all elements in the domain. 

The stationary values of the resulting expression yield the following 
variational statement 


6tt = Z [-fr, P + 4) 6ct) )dfi 

e ■- Q e ^,x ^,x ^,y ^,y 


+ /s f^ 6<pds] = 0. 

2e 


( 20 ) 


We 

element 


then choose bilinear shape 
the interpolation 


0 

functions N^(x,y) 


so that within each 


(f)^(x,y) = N^(x,y)({)^ = ,(i=l,...,m) 


( 21 ) 


0 

holds, where (j)^ are the nodal values of velocity potentials and m is the 

number of element nodal points. Substituting equation (21) into (20), we 
obtain the following system of equations 


K ^ ^ (22) 

where 

K=Z/^p(N N^ + N N '^)dQ (23a) 

- e ^^^^e -,x-,x -,y-,y' 

f = Z / f NdS (23b) 

C ® 

e S 

26 


^ . (23c) 

e — 


Four-node isoparame trie quadrilateral elements shown in figure 2 with 
the bilinear shape functions (ref. 35) 

N® = i (1 + (1 + nn.) (24) 

are employed for this investigation, where ^,r) are the natural coordinates 
and denote the nodal point coordinates of the parent square element in 
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figure 2a, With the same shape functions, transformations in the form 


e 


X 


N?x. 
1 1 


e 

y 


N?y . 

1 1 


( 25 ) 


are used in mapping the parent element to a general quadrilateral shown in 
figure 2b, where x^ and are the global Cartesian coordinates of the nodal 

points for the transformed element. 


A selective Gaussian integration scheme is employed for evaluating the 
element matrices in equation (23a) . The density is evaluated only at 

the centroid of the elements, while a two-point Gaussian quadrature per 
direction is utilized for the area integration of the shape functions and 
their global derivatives. Since the local density in an element is nearly 
constant in the flow direction, no sacrifice in accuracy is made. This 
leads to significant savings in computations because the area integrations 
of global derivatives need be performed only once throughout the entire 
nonlinear iterations. 


Pseudo-Unsteady Formulation 

To obtain a steady-state solution to the nonlinear equation (1), a 
pseudo-time formulation is considered in the following form: 


— (pd) ) + — (Pd) ) + (p(j) ) + (p(|) ) = 0 

CO ^^^,xt ,x to ^^^,yt%y ^^^,x%x ^^^,y ,y 


(26) 


where At is the pseudo-time increment and to is a relaxation factor or a 
damping coefficient. The weak form of the above equation is: 




+ (P^ J ^ ]6(j)(x,y, t)dS^ = 0 

9^ yj yj 


for Integration of equation (27) by parts yields: 

6tt = -/^p(— (j) ^6(j) + — (f) 6(j) +4) 64) 

^2 CO ^,xt ^,x to ,yt ^,y ^,x ^,x 


+ 4) 64) )dn + f f6cj) ds = 0 

,y ,y ^2 


where f = p4) is the mass flux specified on the boundary S^, 
, n ^ 


(27) 


(28) 
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Finite element discretization of equation (27) yields the following 
pseudo-unsteady system of nonlinear equations: 


Ate (j) + K (}) = f 


( 29 ) 


where 



(30) 


is the damping matrix of the above time dependent problem. Here oj is a 
relaxation parameter controlling the stability of the numerical integra- 
tions. Using backward - differencing in time, the following step-by-step 
numerical integration scheme for equation (29) is obtained 

^n+l^n+l ^ ^n+1 ^3 


with 


, n+1 ^n+l , , T , n 
(f) =030 + (1-0)) 0 


(32) 


where n is the pseudo-time step. 


Equation (31) is nonlinear and a "good" estimate of the coefficient 
matrix is necessary at each time step to assure the stability of the nu- 
merical integrations. It has been shown in reference 15 that 

an appropriate procedure for evaluating and consequently at a 
time step n is to use 


n+1 

P. 



when the local flow is subsonic, and use 


(33) 




A II 

- a As p 
e e e 

,s. 


(34) 


when it is supersonic. In the above equation s is the streamline direction. 

As is the element size in the direction of s, and a is known as the co- 
e e 

efficient of artificial viscosity. 

The extrapolation scheme of equation (33) takes into account the 
elliptic nature of the equations in the subsonic region whereas equation 
(34) takes into account the hyperbolic nature of the equations in the 
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supersonic region. For sufficiently streamlined elements and uniform 
grids, equation (34) can further be expressed as follows 


n+1 


CL 

e 


n , 1 \ n 

p + (1-a )p 
eu e e 


( 35 ) 


where p^^ is the mass density of the nearest element at the upstream side 

of e. This corresponds to the well known "upwinding" technique used in 
finite differences (ref. 19). In the following section, the amount of 
artificial viscosity a and the relaxation factor oo to be used will be 

determined according to the stability considerations of the numerical 
integrations. 


Obtaining an accurate and efficient solution of the potential flow 
equations in the transonic regime has been a major concern during the recent 
years. A popular approach has been to extend the finite difference pro- 
cedures originally developed for the solution of subsonic flows by adding 
artificial viscosity in the supersonic region. In reference 19, various 
types of iterative procedures for solving the same equations are reviewed. 
ADI, SOR and multi-grid methods as well as higher order relaxation schemes 
have been applied to improve the convergence properties of transonic flow 
equations (refs. 8,19,20,21,23,26,31,32). The finite element method, on the 
other hand, has also been applied to the solution of elliptic problems by 
many investigators for several years. The extension of this method to 
transonic flows has been accomplished basically by following the same pro- 
cedure, i.e., by adding artificial viscosity to the original system of 
equations (ref. 15). 


The Numerical Integration Technique 
for Pseudo-Time Problem 

The numerical integration scheme based on a direct solution of the 
system defined in equation (31) has been employed with success for the 
solution of transonic flow problems in reference 15. However, since the 
left-hand-side of equation (31) is nonlinear, a full decomposition of the 
coefficient matrix K is needed at each step which makes the scheme com- 
putationally somewhat inefficient. For efficiency considerations the time 
integrations in equation (31) are modified as follows: 

K« + (k" - 

with 

,n+l n+1 , .,n 

£ = + (l“0))j^ 


(36) 


( 37 ) 
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where 


K° = K(p ) and p = p. (38) 

— — ''^ 00 ^ ^00 in 

are constant for all time steps and is evaluated using p^"^^ as given in 

equations (33) and (34) . Hence, the decomposition of K° is needed only for the 
first time step and the subsequent solutions can be obtained with relatively 
inexpensive forward and backward substitutions. We shall refer to the 
scheme in equation (31) as the variable coefficient scheme (VSC) and the one 
defined by equation (36) as the constant coefficient scheme (CCS) . 


Modeling of Shocks 

Considerations for convergence and accuracy in the shock region rely 
strongly on the modeling of the shock. Either shock capturing or shock 
fitting techniques can be employed to model the shock. Shock capturing 
techniques have recently become more popular due to the generality of 
locating the shock and the simplicity of the modeling. However, to attain 
convergence and stability around the shock, the choice of the computational 
grid and the artificial viscosity parameter requires considerable effort 
with either shock capturing or shock fitting methods. The basic consider- 
ations in analyzing the stability and the convergence of the numerical 
scheme in the supersonic flow do not remain valid since the problem in- 
volves a shock discontinuity during the transition from supersonic to sub- 
sonic flow. To localize the smearing that occurs while modeling the sharp 
changes at the shock, the computational grid often needs refinement. More- 
over, to stabilize the high frequency oscillations in the vicinity of the 
shock, artificial viscosity is sometimes added even to portions of the sub- 
sonic flow immediately following the shock (ref. 26). The main disadvantage 
of this approach, hence, is the need for judicious modeling in this area. 

On the other hand, if one attempts to use the shock fitting procedure 
directly, convergence difficulties are reduced considerably. However, one 
is faced with the basic problem of locating the shock accurately. 

In the present formulation, the potential function (p is continuous 
across the element interfaces but the derived quantities p,p and a are 
discontinuous. Supersonic elements are treated by adding artificial vis- 
cosity as discussed previously. While the flow is changing from supersonic 
to subsonic, the continuity of the mass flux in the form 


-+, + ~ 

p(p=p(p 

• n , n 


(39) 


is satisfied as a natural consequence of the variational formulation, where 
+ and - denote the upstream and downstream sides of the interface, and p is 
the modified density defined by equation (34) . Hence, providing that the 
interfaces are sufficiently aligned with the shock lines, the shocks are 
allowed to appear between the elements, and the continuity of the actual 
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mass flux 


(40) 


+ , + — , ■ — 
P 4) „ = P „ 

»n ,n 


is satisfied in the limit, as the artificial viscosity parameter approaches 
to zero. For finite values of equation (40) can be imposed as an addi- 
tional constraint to the variational form, as will be discussed on page 17. 
However, the effect of violating equation, (40) is found to be insignifi- 
cant for small values of a . 

e . ■ ■ 


Since, across the shock, the influence from downstream to upstream is 
still present, the above is basically a shock capturing scheme. The scheme 
requires grid refinement at the shock in order to localize the smearing 
effects due to downstream influence. Moreover, high frequency oscillations 
appearing in the shock vicinity require relatively more artificial viscosity 
in this region to damp out such oscillations. The resulting solution at the 
shock is extremely sensitive to the way. these difficulties are handled, 
consequently, no unique solution is readily available. To remedy this situ- 
ation, we propose to use a shock fitting technique for checking the unique- 
ness of the obtained results. 


In the shock fitting procedure, the shock is also assumed to occur at 
the element interfaces. However, after a shock is detected as a change 
from supersonic to subsonic flow between two adjoining elements, these 
elements are uncoupled, and appropriate boundary conditions are imposed 
on both sides of the shock (ref. 15). Although this technique is computa- 
tionally straight-forward to implement, as in all other shock fitting tech- 
niques, it relies strongly on the prior determination of the shock position. 
For this reason a combined version of shock capturing and shock fitting 
techniques is employed. 

a) Shock Capturing 

In the case of shock capturing, no modifications are made in the basic 
form of finite element equations. Although shocks can still occur as dis- 
continuities between elements, smearing is present due to the coupling be- 
tween the subsonic and supersonic elements neighboring the shock. An addi- 
tional constraint equation is introduced along the shock, since the con- 
servation is no longer guaranteed at the shock interface because of up- 
winding. The variational functional in equation (18) is modified as 
follows : 


= Z TT + AZ (pV - P 4>~ )^ds 

e ^ e^,n ^ 

e f f ’ 


(41) 


where A is a penalty constant, is the common edge of the elements on the 

shock and the symbols (+) and (-) denote each side of the shock. The im- 
plementation of this concept is simply achieved by checking the adjacent 
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elements in subsonic and supersonic flow regions. The element equations 
are then modified with the inclusion of the line integral in equation (41) . 
Details of such interface elements are given in reference 2. 

b) Shock Fitting 

In the shock fitting procedure, the shock is again assumed to occur 
at the element interfaces. However, when a shock is detected as a change 
from supersonic to subsonic flow between two adjoining elements, these 
elements are uncoupled as shown in figure 3. This results in a physical 
separation of subsonic and supersonic regions of flow on the computational 
grid. 



Figure 3. Shock fitting model 

A natural boundary condition is imposed on the shock line of the super- 
sonic element in the following form: 


^n+1 ^n . ^n 

f = f - a As f 
e e e e e , s 


(42) 


where 


, +, +.n 

f = (p 0 ) 

e e ,n 


(43) 
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The value of at the shock points are in turn placed as forced bound- 
ary conditions 


<t. = 4 

at the corresponding nodes of the subsonic element along S^. Conservation 

is once more satisfied through the use of the penalty function technique, 
but this time only for the downstream element, by modifying the functional 
as follows: 


A 

7T 

e 




(45) 


In equation (45) no variations are taken with respect to (+) terms since 
they are assumed to be known from the upstream conditions. No additional 
node points are created. This requires the modification of the element 
equations along the shocks. In this way, fundamentally different charac- 
teristics of subsonic and supersonic flows are separately modeled. 


Since, as it will be shown that the scheme is uniformly convergent 
in a wider range of initial conditions for higher artificial viscosities, 
shock capturing solutions are first obtained by selecting a large arti- 
ficial viscosity parameter a^. Following this, a series of shock capturing 

solutions are obtained by gradually decreasing until convergence diffi- 
culties are detected. Using the previously converged shock capturing 
solutions as an initial guess, an attempt at shock fitting is subsequently 
made to check the shock position. If the initial shock position is correct 
and the artificial viscosity content is optimum, the influence from down- 
stream to upstream at the shock is minimal in the shock capturing. Hence, 
both techniques yield essentially the same solution beyond which no further 
decrease in the artificial viscosity is possible. 

If the initial position of shock is incorrect, the shock moves to a 
different position after fitting, generally towards the downstream direction. 
The developed fitting scheme differs from the classical fitting procedures 
for which the position of the shock remains fixed. In the present case, 
however, when the conditions in equations (42) and (44) are applied on a 
particular flow configuration, a new shock position is free to develop after 
each iteration. The computer code automatically detects the new shock 
position and applies fitting once more. In that case the shock capturing 
is again tried by lowering the artificial viscosity, and the above process 
is repeated until a unique solution is obtained. 


A novel feature of the present formulation is the ease with which 
a combined version of both techniques is utilized. At any point 
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•during the iterations, the solution procedure can be switched from one 
scheme to another to improve efficiency and accuracy without altering the 
computational grid. It should also be mentioned that the adaptation of the 
grid to the direction of an oblique shock which has been successfully ap- 
plied by using shock fitting techniques can also be implemented in the 
developed procedure. 


Computational Grid, Periodic Boundaries 
and Kutta Condition 

A typical computational grid shown in figure 4 is employed here in 
modeling highly staggered cascades. The broken lines on the grid denote 
the far field inlet and exit boundaries of the solution domain shown in 
figure lb. The extensions at the inlet and exit boundaries are used to 
maintain a rectangular outer domain instead of a parallelogram for effi- 
ciency and ease in the grid generation. Uniform farfield flow conditions 
are imposed in these extended regions. The use of rectangle-like elements 
in the grid has the advantage of providing grid lines which are aligned 
better with the direction of the flow as well as the direction of the 
shock lines. Moreover, with this grid severe element distortions can 
readily be avoided at the rounded leading and trailing edges. 

Along the periodic boundaries, every node need not have a matching node 

on the opposite periodic side. The matching periodic nodes are identified 

with solid circles in figure 4 for which the periodicity constraint given 

by equation (15) is imposed. Equation (15) is satisfied by assigning the 

same equation number to each pair of the periodically located nodes and 

modifying the right-hand side of the element equations on the boundary line 

B-D in figure lb, with the potential difference q. (h CosS, + h SinS . ). 

in^ X in y in 

The nodal values of the intermediate periodic nodes are computed by inter- 
polations from the values of the nearest pair of nodes. The final peri- 
odicity is then attained iteratively. This scheme allows a flexibility in 
placing the periodic nodes, which leads to the generation of more efficient 
grids with refinements only in the regions of interest. 

A symmetric frontal solver (ref. 22) is employed in solving the system 
of equations. In this solution procedure, only those equations that are 
currently required for the elimination of a specific nodal point variable 
need by assembled and saved in the high-speed storage. As opposed to band- 
ed solvers, the effectiveness of the scheme does not depend on the nodal 
point numbering sequence. Hence, the periodic boundary conditions can be 
imposed conveniently by assigning the same equation numbers to each pair 
of periodically located nodes. This may be accomplished without renumbering 
the grid points, and also with no loss of computational efficiency. 

If the exit flow conditions are not known a priori, the Kutta condition 
at the trailing edge is iteratively satisfied through the following proce- 
dure. The circulation around the airfoil is determined from a previous 


18 



iteration, using the computed velocity potential values of the two un- 
coupled nodes at the sharp trailing edge, as follows; 

X I ^ (46) 

where + and - denote the two sides of the slit which extends from the 
trailing edge to the exit boundary. This value of circulation is then 
imposed all along the slit as a constant discontinuity in velocity poten- 
tial for the new iteration. The procedure is repeated until an overall 
convergence is reached and the Kutta condition 

q I = q I (47) 

at the trailing edge is satisfied. 
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CONVERGENCE CHARACTERISTICS OF THE EQUATIONS 


Convergence of Variable Coefficient Scheme (VCS) 


The computational procedure outlined thus far involves appropriate 
choice of in equation (34) and co in equation (37) in obtaining accurate 

and convergent results. To investigate the convergence characteristics of 
the equation, to a steady state solution, a simple one- dimensional model 
will be used, and since the convergence limits obtained from a single 
element model provides an upper bound to the general system, only one 
element will be considered. Hence, for an element e and for y = 2, 
equation (31) takes the form 


ns2 


{K^ - (q:) 

0 


+ a Ax [(q’^)^l }N 
e e e , x — . 


,,T ^n+1 = 0 
N * 

X— ,x— e 


(48) 


If one assumes that the estimate to the nth solution is close to the 

steady state solution 0 or q in element e so that 

^ e 




ll,M » I 


(49) 


~n , ,~n 

q = q + Aq 
e e e 


q^l >> l^q^ 


where 


denotes a suitable vector norm and 


(50) 

denotes an absolute 


value, then equation (48) becomes 

^2 ^2 , A , /„2 


[K^ - q^ + Ct Ax (q^) ]N 

e e e e ,x — , x— , x ^ 


= [2q Aq’^ - 2a Ax (q Aq“) ]N N 

e e e e e e ,x — , x— , x^ 


n. 


Recognizing that 


T .;^n+l .~n+l j A 


(51) 


(52) 


equation (51) can be rewritten as follows: 


[K - q + a Ax (q ) ]N Aq 
*■ e e e e ,x — ,x e 


n+1 


= [2q^ - a Ax (q^) ]N Aq^ - 2a Ax q^N Aq^ 

■■ ^e e e e ,x-^— , x e e e e— ,x e,x 


(53) 
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By using equation (26) it may further be shown that 


. - Aq^ 2q^ - 2a^ - 2a Ax (q^) 

At ^e ^e N ; ■ e e e e ^e ,x * n„ 

— 7— »x = ^ Aq N 

^ o2, a/-2n e— , X 

2a + a Ax (q ) V 

e e e e ,x 


2a Ax q^ 
e e^e 


2a^ + a Ax (q^) 
e e e e ,x 


Aq^ N 
e,x-,x 


( 54 ) 


where 


a" = (K^ - q^)/2 
e e 


(55) 


The above difference relationship implies that the error Aq (x,t) 
satisfies the following differential equation: 


* n , . n, ^ . ri 

Aq + cAq = bAq 

e , t e , X ^e 


(56) 


where for small amounts of artificial viscosities 


. , .2 

c = a Ax M -r— , 

e e e At * 


b = (M^ - , 


(57) 

(58) 


and 


is the local Mach number. This model error equation will be employed to 
describe the basic characteristics of the pseudo-time integration procedure 
in the different flow regions. 

a) Convergence in Subsonic Flow, c = 0 : 

In the case of subsonic flows, with no added artificial viscosity, the 
solution of this problem requires accurate and efficient integration of the 
reduced system 


Aq^ = bAq^ , b < 0 
e , t e 


(59) 
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When the local flow is subsonic, i.e., when M < 1, even with zero artificial 

viscosity a , the scheme converges to a steady state provided that co > 0. 

0 

The solution of equation (56) in this case is 

= Aq^exp[bAt] (60) 

which represents the response of a first order damped system. The rate of 
convergence depends on both the relaxation factor iii and the local Mach 
number. In addition, to assure the stability of the numerical integration 
procedure in time, the bounds for O) become 

0 < 03 < 2/(1 - M^) (61) 

Equation (61) shows that for compressible flows, the over-relaxation can 
be higher as a function of local compressibility compared to incompressible 
'flows. Multi-grid methods or the use of over-relaxation schemes mainly aim 
at increasing the convergence of the iterations for small values of b, 
which is typical of elliptic systems. 

b ) Convergence in Supersonic Flow, c 4 Q • 

When the local flow is supersonic and no artificial viscosity is used 
(i.e., = 0) the solution of equation (56) becomes divergent as implied 

by equation (60). However, when a positive artificial viscosity is intro- 
duced, the characteristics of the system change and it becomes a hyperbolic 
system. The solution of equation (56) for > 0 can be written as 

Aq^"*”^ = Aq^ (x-cAt) exp [bAt] (62) 

This solution indicates that since the errors in this case cannot be 
damped, they are convected in the flow direction with a wave speed c 
which is linearly proportional to the amount of added artificial viscosity 
and the relaxation factor o). While propogating however, the errors are 

magnified since b > 0. The rate of magnification depends on the values 

and 03 as defined in equation (58) . 

By changing the behavior of the equations into a convective form with 
the addition of artificial viscosity, convergence to a steady state is 
guaranteed in equation (62) for a fixed point in the supersonic pocket. At 
a particular time step during the iterations, typical variation of errors 
in the supersonic pocket is as shown in figure 5. Errors at any time step 
have increasing magnitudes in the flow direction and are largest around the 
shock. As the flow reaches a steady state at the subsonic upstream, Aq^, 

at the subsonic point converges to zero. The convergence at an element 
depends on the variation of Aq^ at the supersonic element until all the 
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errors are convected through. It will be shown subsequently that this will 
depend on how far downstream the element is located and on the number of 
elements placed in the supersonic pocket up to that point. 

The above analysis of the convection of the errors in the supersonic 
pocket is general for an "artificial compressibility" scheme, whether it is 
based on a finite difference, finite volume or a finite element method. The 
present study provides a rigorous analysis of the problem. The solution 
scheme for transonic flows should be suited to account for the propogation 
of the error in the flow direction. In the application of the present 
scheme, it was observed that the convection of the errors is the deter- 
mining factor in the convergence rate of the solutions. To illustrate this, 
a sample case is tested for a cascade problem (NACA 0012 airfoils, with 0 

stagger angle, h/C = 3.6, 3- =0, M. = 0.8). The variation of errors for 

07 in in 

each time step are plotted in figure 6 for all of the supersonic elements on 
the airfoil surface, following the numerical integration scheme in equation 
(31) for different values of and oo. In these figures m^ denotes the 

element location at the supersonic pocket; m^ = 1 being the first supersonic 

element past the sonic line. As can be seen, when the flow is subsonic for 
an element, solutions exhibit uniform convergence described by equation (60). 
As soon as the element becomes supersonic, an artificial viscosity is added 
and the errors are convected with a wave speed linearly proportional to the 
artificial viscosity from equation (57) . Also from the same equation, the 
change in the relaxation parameter, 03 , increases both the damping factor, 
which is negative in this case, and the wave speed. Therefore, the errors 
are magnified as they are convected along the flow direction until they reach 
the shock. The convergence of the solution is complete when all errors gen- 
erated by the upstream elements have been convected. This makes the shock 
region critical for convergence and determines the stability characteristics 
of the scheme. 


Convergence of Constant Coefficient Scheme (CCS) 


Following the same procedure presented in this chapter, we consider a 
one-dimensional version of equation (36). For a single element e and y = 2, 
equation (36) takes the form 


{K^ - q^} N N 


(q^ + (q^)2 - a Ax [(q"")^] } N N = 0 . 


(63) 
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Figure 6. Convergence to the steady state in the 
supersonic pocket (Continued) 





If the estimate to the nth solution is close to the steady state solution 
(p in element e, i.e., if 


= 0 + A$^ with I |(f) II » j IA^’^I 


--e 


(64) 


and 


= q + Aq^ V7ith q >>lAq^l , 
e e e e ' n ' 


(65) 


then equation (63) becomes 
■2 ^2 


[K^ - q^]N N 
*■ , X— , X — e 


= [-2q^ + q^ - a Ax (q^) ]N N "^Acj)^ 
oo e e e e ,x — , x—,x ^ 


+ [2q Aq’^ - 2q Aq^ - 2a Ax q Aq’^jN N (}) 

'■ ^e e “ oo e e e e — ,x— ,x-^ 


( 66 ) 


Recognizing that 


T.^n+l .~n+l j T, 

N A©^ = Aq and N A = q , etc., 

— , X ^ e — ,x-^ e 


(67) 


equation (66) is simplified as follows 


[K^ - = 

= [-2q^ + 3q^ - 2a Ax (q^) ]N Aq^ - 2a Ax q^N Aq^ 

>- ^oo e e e ,x — ,x e e e e— ,x e,x 


( 68 ) 


Using equation (37) and the identities: 


2a^ = - q^ ; 2a^ = - q^ 

e e °° 


(69) 


we obtain 


. n+1 . n 

At a^ *^e ^ N 

(jO ” At 


= [q^ - a^ - a Ax (q^) ]Aq^ - a Ax q^Aq^ N 
^ e e e e e ,x-* e— , x e e e ^e,x— ,x 


(70) 
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The above difference relationship implies that the error again satisfies 
the following differential equation 


Aq"^t + 


(71) 


where for small amounts of artificial viscosities 

a^ 

c = a Ax ^ ^ 

e e e a 2 At 


(72) 


b = 


(m" - 1) ^ 
a 2 ^ At 


(73) 


Comparison of the above stability analysis with the results obtained 
for VCS in equations (56-58) shows that in this case the coefficients b and c 
are both magnified with a factor of a^ will be shown in the follow- 

ing section that the choice of artificial viscosity, from equations 

(57,58) or from (72,73) provides the same boimds for convergence. It can 
be concluded that both CCS and VCS exhibit similar convergence character- 
istics, but CCS is computationally more efficient. 


Choice of the Artificial Viscosity Distribution 

One of the most important factors in determining the efficiency of the 
solution procedure has been the choice of the artificial viscosity distri- 
bution. The finite difference schemes have traditionally employed a fixed 
artificial viscosity distribution for a particular grid. For each grid 
point, this distribution is determined only as a function of the local Mach 
number. Thus, while for rough grids the convergence is fast, for finer 
grids it slows down due to the reduction in the artificial viscosity. 

The basic steps of the finite element scheme developed here are as 
follows: i) choose a relatively high artificial viscosity, ii) obtain 

convergence to a steady-state, iii) reduce the artificial viscosity, 
iv) repeat the same steps until a solution with sufficiently small arti- 
ficial viscosity is obtained. This procedure has proven to be very re- 
liable; one can always converge to a stable solution which is also checked 
for accuracy by the shock fitting technique. It is similar to the appli- 
cation of successive mesh refinement techniques where one is increasing the 
number of mesh points for reducing the artificial viscosity. 
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For a better understanding of the choice of the artificial viscosity 
distribution, equation (56) is further studied. The solution of this 
equation can be written in the following form: 

Aq^ = Aq° (x-ct) exp(bt) (74) 

If an initial solution, Aq°, is assumed in the form of a ramp function, 

Aq° = 0 for x < 0 and x > L 

e 

Aq° = q x/L 0 X L (75) 

e max 

where L is the length of the supersonic pocket, the condition for uniform 

convergence, Aq^ < 0, provides 
e , t 


c 


bx 

1-bt 


(76) 


The condition for uniform convergence requires the following lower limit 
for the artificial viscosity: 


m 


a 

e 


(1 - 


) 


(77) 


l-o)n(l - 






where is the element number indicating its position in the flow direction 

starting from the shock, and n is the number of time step (t = nAt) . As can 
be seen in this formulation, the most undesirable case occurs at initial 
time step (n = 0) for the element neighboring the shock, m^ = N, where the 

length of the supersonic pocket is L = NAx. For both CCS and VCS schemes, 
if only the initial time step is considered, then, equation (77) reduces to. 


a ^ m (1 - — 
e e 

e 


) 


(78) 


In references 15 and 4, a limit for the artificial viscosity was defined as, 

1 „ 1 


a :>^ , 
e 1-r 


(1 - 


-) 




where r is the ratio of error for two neighboring elements, i.e 

r = -7 < 1 . 

e Aq 


(79) 


( 80 ) 
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It can be shown that, for the case of a ramp function as an initial 
solution, equation (79) reduces to equation (78). Also, if the propa- 
gation of the errors to the neighboring elements is neglected, for m = 1, 
equation (78) becomes, 

a I - 1/m’ (81) 

e e 


which is the commonly employed expression. 

The above solution of the model equation indicates several important 
considerations in the selection of an artificial viscosity distribution: 

. Use of higher values of is necessary at the initial 

stages of the numerical integration procedure. However, 
lower values of artificial viscosity are permissible at 
later stages as can be seen from equation (77). 

. Equation (77) also suggests that the amount of artificial 
viscosity should also be determined on the basis of the 
position of an element in the supersonic pocket. 

. As can be seen from equation (79), the amount of artificial 
viscosity necessary for uniform convergence depends on the 
local Mach number as discussed by several investigators, 
e.g. (refs. 19,21,25). 

. In addition, however, the required artificial viscosity 
depends highly on the distribution of the errors in the 
initial solution which is implied by r^ in equation (79) . 

For larger values of artificial viscosity, the scheme is 
uniformly convergent for a wider range of initial conditions. 

. Also, as can be seen from equation (80), when a grid is 
refined and the element size gets smaller, approaches 
unity and the value of artificial viscosity necessary 
for uniform convergence increases. 

Customarily, the non-dimensional artificial viscosity coefficient a 
as defined in equation (81) is employed in the literature as the ^ 

measure of the amount of artificial viscosity. However, for a comparison 
of accuracy of solutions obtained with different computational grids, one 
should use equal values of 

a = a As (82) 

e e e 

rather than for both grids as a measure of the total artificial 
viscosity employed in equation (34) where As^ is the size of the element 
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employed in equation (34). If one compares solutions obtained from 
different grids by using same values of a^, the results will be different 

not only due to the discretization errors but, also, due to the different 

quantities of artificial viscosity a^. In fact, the latter becomes the 

main source of error in the supersonic region. For the same reason, in 

the case of non-uniform grids a correction is needed on a to maintain a 

e 

uniform artificial viscosity distribution. This correction is achieved 
simply by using a modified artificial viscosity coefficient as follows: 


a' 

e 


a 

e 


As 


av 


As 

e 


(83) 


where is an average element size in the supersonic pocket, and As^ is 

the average of the lengths of elements e and eu; eu indicating the element 
number at the upstream of element e. Consequently the total artificial 
viscosity used in each element becomes 

As As 

a = a As = a (84) 

e e e e - 7 — 

As 

e 


Thus a uniform distribution of a rather than a is obtained for the 

e e 

irregular grid. 


a) Von Neumann Analysis of Convergence : 


For a better understanding of the convergence characteristics of the 
equations in the supersonic region a Von Neumann type stability analysis 
of the error— equation (56) is performed after expressing a difference 
equation in time as follows: 


. n+1 . . r A.- cAtv A n. cAt . n 

Aqe = (A + bAt - Aq^ + ^ Aq^^ 


(85) 


Assuming that a supersonic region of length L is divided into equally 
spaced elements of size Ax^ = L/N, the Fourier components of the error at 

a mesh point k and a time step n can be written in the form: 


where 


Aq’^ = exp(le0, ) (86) 

Ax 

I = = TTk , k = 1,2 , N+1, 

ttAx 

min0^ = ^ ^ , max0^ = it (87) 
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Substitution of equation (86) into equation (85) yields the following ampli 
fication factor for the magnitude of error: 


V 

V 


n+1 

n 


1 + bAt - 


cAt cAt 

Ax Ax 
e e 


(CosGj^ 


ISinGj^) . 


( 88 ) 


The polar diagram of the amplification factor is shown in figure 7. 



Real 


Figure 7. Polar diagram 

The two necessary conditions for uniform convergence as observed from 
figure 7 are: 

1) R = <1 or a GO < (89) 

e M 

e 

2) S = 1 + bAt - <1 or, a > 1 ^ (90) 

^ 2 . 

e M 

e 

Equation (90) provides a restriction on the relaxation parameter. At 
high artificial viscosities, as the propogation speed of errors are high, 
the relaxation parameter must be decreased for retaining stability of the 
numerical integrations. Since in subsonic flows 

R + S = 1 + bAt < 1 (91) 
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the scheme is always convergent. For supersonic flows, however, the 
inequality in equation (91) is not satisfied. Therefore, the scheme is 
not uniformly convergent unless 


0 . 

min 


> 0 


cr 


( 92 ) 



1 

1 + 0)(M ^ - 1) - M ^ 
e e e 


(93) 


The above equation indicates that the convergence at a particular super- 
sonic point depends on the number of elements, N^, counted from an upstream 

point at which an error starts propogating. As the upstream points converge 
to a steady state, for a particular point decreases requiring less arti- 
ficial viscosity. At a certain time during the iterations the convergence 
characteristics change along the downstream direction based on the upstream 
conditions . 


This analysis demonstrates again the importance of accounting for the 
convective nature of the errors in the design of both the solution proce- 
dures and the computational grids. 


ACCURACY OF THE SOLUTIONS 


Introduction 

The accuracy of the transonic flow solutions depends on the correct 
modeling of the several flow regions with different properties. In the 
analysis of flow through a cascade, basic regions of interest are the 
supersonic pockets, the shock lines, and the leading and trailing edges of 
the airfoils. The application of shock capturing techniques, with added 
artificial viscosity for obtaining stable solutions, at the same time 
complicates the accuracy problem. One can introduce "wiggles" or "exces- 
sive smearing" to the solution through the artificial viscosity added to 
the system. It becomes quite difficult to analyze the accuracy of the 
obtained solution as a function of the number of node points if the amount 
of artificial viscosity changes with the grid refinement. The accuracy 
problems can be categorized in the following three main groups: 

(a) Discretization errors due to high gradients in pressure distribution. 
These errors are important generally around the leading and trailing edges 
of the airfoils where the flow is subsonic and the equations are elliptic. 

(b) Errors due to the artificial viscosity distribution. Since the con- 
verged solution at the supersonic pocket is obtained by a specified arti- 
ficial viscosity distribution, the difference from the inviscid flow can 
be considered as an error. The minimization of this type of error requires 
proper application of the artificial viscosity distribution as a function 
of local flow and the computational grid. 

(c) Errors in modeling shock discontinuities across which the flow 
equations change their characteristics. 


Computational Grid for Subsonic 
Flow Regions 

In the case of subsonic flows, equation (1) is elliptic and the 
finite element formulation involves the calculation of the local maxima 
of the variational functional in equation (18) . The convergence to a 
steady state is guaranteed from the integration of equation (26) as shown 
in equation (60) . The rate of convergence on the other hand can be im- 
proved by using relaxation as implied by equations (60) and (61) where 
relaxation factor w is determined as a function of the local Mach number. 

In terms of designing an appropriate grid for subsonic flows, one has to 
be mainly concerned with the accuracy and the faster convergence of 
iterations. It has been shown during recent years that the computational 
grids can be important in improving the convergence of the relaxation tech- 
niques expecially, for elliptic equations (ref. 5). 
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From the mathematical standpoint, the finite element discretization 
of equation (1) requires the calculation of the pressure integral over 
each element area individually. The discretization errors for the entire 
flow becomes the sum of the discretization errors in each element. For a 
given number of elements, the grid has to be designed to minimize the total 
discretization error. It sould be noted that, because of the elliptic 
nature of the equations in the subsonic region, local discretization errors, 
due to boundary geometry generally remain local and can be improved by a 
gradual refinement of the mesh in the immediate problem region. The 
accuracy problem related to choosing computational grids has also been 
studied by finite element researchers for several elliptic problems, e.g., 
(refs. 33 and 34). 

In designing a computational grid for subsonic flows, one has to refine 
the grid at regions of high pressure gradients, such as leading and trailing 
edges of an airfoil, until the results reach a desired accuracy level. 

Local modif ications of the grid have generally localized effects on the flow 
configuration. When the grid is refined, the convergence follows the char- 
acteristics of the model error-equation in equation (59) . 

To illustrate some of the points discussed, incompressible flow for a 
Gostelow airfoil cascade was analyzed. The simple sheared grid shown in 
figure 8 was used. The results shown in figure 9 are representative of the 
problems encountered in the solution of subsonic flows in cascades. At the 
leading edge of the airfoil a near singularity occurs. In this area of 
high pressure gradients local discretization errors dominate the solution. 
This is attributed to the high distortion of elements around the leading 
edge as shown in figure 10a. Figure 10b illustrates the simple adjustments 
made on the initial grid to improve the situation. As can be seen from 
figure 9, the improvements made by the change in the grid remain local. 


Designing Grids for Transonic Flows 

The design of computational grid in transonic flows requires several 
considerations. The accurate representation of the subsonic leading edge 
becomes more important in predicting the supersonic pocket as the sonic 
line develops closer to the leading edge. The accuracy of the solution 
is most critical at this region. As discussed previously for subsonic 
flows, an appropriately refined grid is necessary for the modeling of high 
pressure gradients. Similarly, accurate modeling of the trailing edge be- 
comes important since it affects the accuracy of the exit flow conditions. 

The source of difficulties for modeling the supersonic pocket in a 
cascade is quite different. As it will be shown with numerical examples, 
the problem is more that of convergence than of accuracy. The classical 
concept of refining a grid for better accuracy has to be very carefully 
applied in this case. The choice of artificial viscosity and relaxation 
parameter as well as the modeling of the shock are also strongly dependent 
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Figure 8. Computational grid for Gostelow airfoil 





Figure 9. Comparison of results - 
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Figure 10. Grid adjustment at the leading edge 
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on the computational grid. When one refines a computational grid for 
accuracy, the efficiency of the solution requires a new set of these 
parameters . 

To illustrate some important aspects in analyzing the transonic flows, 

a cascade of NACA 0012 airfoils with a pitch to chord ratio of 3.6, and the 

flow conditions M. = 0.80, 3- = 0.0°, 3 “ 0.0°, and (x = 0.0° was con— 

in in ex s 

sidered. Four different grids shown in figure 11 were used for this non- 
lifting case, where, due to symmetry in the cascade geometry only one half 
of each grid is drawn. There are 10 elements over the airfoil in the 
coarse grid A, 20 elements in the intermediate grid B, and 28 elements in 
the fine grids C and D. Grid A and B have 12 elements between each airfoil, 
whereas, grids C and D have 28 elements. There are a total of 312 elements 
in grid A, 480 elements in grid B, and 1152 elements in grids C and D. 


Figure 12 shows the variation of pressure coefficients C over the 


chord length C for the coarse grid A. A series 
were obtained by employing different amounts of 
this grid. Artificial viscosity was added only 
as given in equation (34). In order to account 
gation of errors to the neighboring elements, a 

is employed throughout the supersonic region so 
in equation (79) is calculated as 


of solutions presented 
artificial viscosity for 
to the supersonic elements 
for the effect of propo- 
nondimensional constant y 

e 

that the artificial viscosity 


a 

e 



e 


(94) 


Furthermore, for comparing the results obtained from different grids the 
artificial viscosity multiplier is scaled such that 


y 


e 


As 


As. 


(95) 


is used in measuring the artificial viscosity content for each case, 
where As^ is the average mesh size in the supersonic pocket of grid A. 

A series of solutions were sequentially obtained by using the shock 

capturing technique for y^ = 3.0, 2.0, 1.5 and 1.0. The shock developed 

between the same elements in all cases although its position was not correct. 
An overshooting for low values of the artificial viscosity also indicates 
this inaccuracy in the shock location. A shock fitting was then applied to 

jSt * 

the solution at y^ = 1.0 and subsequently a solution at y^ = 0.8 was obtained. 
In figure 12 only the results of y^ = 0.8 solution is presented since they 
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Figure 12. Comparison of results of Grid A (Coarse grid) - Capturing and Fitting 




are close. The position of the shock moved downstream after the shock 
fitting by one element and the overshooting disappeared. When the fitting 
scheme was switched back to capturing, the solution did not change. 

The observed behavior in these solutions verifies the convergence 
characteristics of the equations as discussed previously in this report. 
During the integration of equation (29) it was observed that the converg- 
ence was obtained in the flow direction in a sweeping manner. The im- 
portance of the path of convergence is evident in these results since the 
location of the shock is determined initially .for a high artificial vis- 
cosity distribution. For small artificial viscosities, as in the case of 
* 

u = 1.0, several solutions can be obtained depending of the shock 
e 

location. The overshoot at the shock region can be explained from the 
inspection of equations (79-80). As r converges to unity, the amount of 
artificial viscosity required is considerable magnified. In the case of 
shock fitting, the uncoupling of the flow regions and the introduction of 
the boundary conditions in equations (42-44) eliminates the overshoot and 
moves the shock in the downstream direction. 

Also shown in figure 12 is the effect of satisfying the mass conser- 
vation across the element interfaces at the shock as defined in equations 
(41) and (45) . It is observed that for high values of artificial viscosity, 
introduction of a constraint for mass conservation modifies the singularity 
at the shock tail considerably. However, as the artificial viscosity con- 
verges to zero, the effect is negligible as it is also evident from 
equations (42) and (45) . 

Figure 13 shows a different set of results obtained again for grid A. 

After a shock capturing solution followed by a shock fitting at = 3.0, 

solutions were successively obtained with different values of artificial 
viscosity by employing shock fitting only. Although solutions with very 
small values of artificial viscosity were obtained rather efficiently, the 
position of the shock did not change. 

Solutions with finer grids B and C exhibited similar behavior in terms 
of convergence and will not be repeated any further. Final solutions ob- 
tained from different grids are summarized in figure 14 and compared with 
the finite area solution of Dulikravich and Caughey (ref. 10). It is to be 
noted that the finite element results presented in figure 14 are for ap- 

■k 

proximately equivalent values of y^ in each grid. In comparing different 

A 

grids, y^ rather than y^ should be used as a measure of the artificial 
viscosity . 

As it is observed from results in figure 14, the accuracy for the flow 
at the supersonic pocket is extremely high even for rough grids, although 
the modeling of the shock and its location needs special attention. On the 
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other hand, obtaining solutions with fine grids requires relatively exten- 
sive computational effort not only in terms of the number of degrees of 
freedom but also in terms of the sensitivity of the solutions to pertur- 
bations . 

Results obtained from grids A and B at a relatively high artificial 
viscosity content are given in figure 15 for comparison. It is interesting 
to note that both grids yield essentially the same results. 

Figure 16 shows the results obtained from a nonuniform grid. The 
artificial viscosity in each element was determined based on the element 
size correction defined in equation (83) .. Also, shown in figure 16 is a 
case in which only shock capturing was used with a slightly higher arti- 
ficial viscosity. Excessive smearing near the shock is of interest in this 
nonuniform case. 

Figure 17 shows the case of a lifting cascade with h/C =' 3.6, stagger 

angle of a = 1.0° and the flow conditions M, = 0.78, 3. = 0.0°, 3 = 1.0°. 

s in in ’ ex 

The results obtained using the same arrangement as in grid B are in a good 

agreement with the finite area results (ref. 11). 


Modeling of the Supersonic Pocket and the Shock 

The accuracy of the solution in the supersonic pocket is determined 
in terms of two basic parameters; the artificial viscosity distribution 
and the modeling of the shock. The importance of introducing correct 
distribution of artificial viscosities for irregular grids was discussed 
in reference 15. It was shown that the incorrect artificial viscosity 
distribution may produce local overshoots or smearing in the solutions. 

In all the previous test cases., it was also observed that the sensitivity 
of the solution increased near the shock. For a 45° staggered cascade of 

NACA 0012 airfoils with h/C = 3.6, and the flow conditions M. = 0.80, 

in 

3^^ = 0.0°, 3g^ = 0.0° (Test Case II , page 76) a numerical study was de- 
signed to show this behavior by shock capturing technique. The artificial 
viscosity parameter is defined, following equation (78), in the form: 


a 


e 


u m 
e e 


(1 



(96) 
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Comparing equations (94) and (96) it is seen that 


y 


e 


LI m 
e e 


( 97 ) 


which is a function of the element position in the supersonic pocket, rather 

than a constant as previously employed. In this case is a constant, and 

the three sets of results in figure 18 include = 0.3, 0.4, and 0.6. All 

of these results are stable and the deviations from each other are higher 
at points closer to the shock. The existence of near solutions for differ- 
ent artificial viscosities in shock capturing schemes is apparent^ especially 
around the shock region. 

The results shown in figure 18 illustrate the typical accuracy problem 
observed in the modeling of the supersonic pockets. It is apparent that in 
all cases the smearing across the shocks are negligible. The two elements 
on each side of the shock (one with and another without the upwinding of 
density) provide an accurate modeling of the shock without smearing. On 
the other hand, the artificial viscosity introduced in the supersonic 
pocket provides smearing for the entire supersonic pocket. Thus, the 
strength of the shock is modified due to the artificial dissipation along 
the entire supersonic pocket. These results suggest that one should be 
concerned with the modeling of the entire supersonic pocket rather than 
just the shock for obtaining accurate results. Figures 14 and 18, also 
suggest that even with rough finite element grids, shocks can be modeled 
accurately with either capturing or fitting schemes. 

One final point to be considered is the use of the artificial vis- 
cosity in the subsonic regions past the shock. This is generally referred 

in the literature as the "cut-off" Mach number, (M < 1) (ref. 26). The 

c 

artificial viscosity parameter, a^, then, is defined with the following 
expression : 


a = y Max [0, 1 ] . (98) 

^ ® M^ 

e 

Depending on the strength of the shock, this usually amounts to introducing 

additional artificial viscosity at the shock region. Generally it improves 

the convergence characteristics of the solution. On the other hand, it may 

produce excessive smearing of the shock region. The comparison of two 

solutions obtained with M =1 and M =0.95 values are shown in figure 19 

c c 

for a sample problem (Test Case IIIc, page 82). These results illustrate 
the sensitivity of the solutions to the introduction of artificial vis- 
cosity at the shock region. Considerable smearing occurs in the shock 
region which distorts the pressure distribution. 
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It can be concluded that the convergence rate should be controlled by 
varying the overall viscosity distribution in the supersonic pocket rather 
than introducing additional artificial viscosity around the shock region. 
The local smearing produced by the above described technique provides con- 
fusing results. It should again be noted that the shock fitting technique 
proved to be valuable in testing the accuracy of the solution in the shock 
region. 


Modeling of the Trailing and Leading Edges 

A common feature of the flow equations in these regions is that they 
are generally elliptic. However, if the supersonic pocket comes closer to 
either the trailing or the leading edge, one has to consider the accuracy 
problem there in more detail. For the case of a near singularity in sub- 
sonic flows, the method of solution is well defined. One usually designs 
a radial grid, with degree of refinement varying as a function of the 
distance from the point of singularity. In the case of elliptic problems, 
the effect of this singularity decays in the radial direction. For the 
test cases of NACA 0012 airfoils with blunt leading edges, no particular 
accuracy problems around the leading edges were observed. In order to ob- 
tain sufficient accuracy at the leading edge, only that particular area had 
to be refined. For the test cases, a coarse model of the leading edge was 
sufficient to provide an accurate representation of the supersonic pocket. 

In the case of the trailing edge of NACA 0012 airfoil, however, the 
sharpness of the trailing edge further complicates the problem. Inaccurate 
modeling of the trailing edge strongly affects the magnitude of the pro- 
duced circulation and changes the overall flow configuration. This be- 
havior was especially found to be important for cascades with high stagger 
angles. Figures (20) and (2l) illustrate the sensitivity of the solution 
to the changes in the Kutta condition. When the exit flow angles were 
changed slightly from that predicted by the Kutta condition, for a sample 
problem (Test Case Illb, page 82), considerably different sets of solutions 
were obtained. These examples indicate the necessity of accurately pre- 
dicting the trailing edge conditions for calculating the transonic flows 
in highly staggered cascaded. 
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a) upper surface 


b) lower surface 


Figure 20. Sensitivity of the flow to changes in the exit angle 

(A: 3 = 4.7° - Kutta condition, B:3 =4.3°, = 5.0 ) 





CASE STUDIES 


To illustrate the general applicability of the developed computational 
procedure and to study the importance of the computational grids for im- 
proving accuracy and convergence, several sample cases were tested. The 
first test case, for transonic flows around a cylinder, shows the sensi- 
tivity of flow solutions to computational grids. The second test case 
illustrates the sensitivity of the transonic flow in the cascade in terms 
of position of the shock and the trailing edge conditions. The third and 
fourth cases are the studies of choked and near-choked flow conditions in 
cascades with high solidity. The final test case is for a highly cambered 
airfoil for which flow angles change considerably. 


Test Case I - Isolated Circular Cylinder 

In order to test the applicability of the developed scheme to external 
flow problems, transonic flow around an isolated cylinder is analyzed for 
the following farfield conditions of 

a) Case l4: M, =0.51 

1 in 

b) Case Ib : M. =0.45 

in 

and compared with results of others (refs. 6 and 19). This problem was 
chosen to investigate both the importance of the computational grid in 
blunt bodies, and the effect of imposing farfield boundary conditions at 
finite distances from the body. It has been observed that the problem is 
particularly sensitive to the choice of finite grids and the type of applied 
boundary conditions at finite distances. 

Three different radial grids A, B, and C, employed for this purpose 
are shown, respectively, in figures 22-24. Due to symmetry of the flow, 
only half of the solution domain is modeled in each case. The farfield 
boundaries extend radially to an outer radius at = lOR from the center 

of the cylinder in each case where R is the radius of the cylinder. It is 
determined that the extension of the grid beyond lOR does not affect the 
results, but as will be shown subsequently, smaller regions do change the 
results considerable. 

a) Case la: = 0.51. -Two types of boundary conditions are imposed at 

the farfield boundaries to simulate the uniform flow field for this case. 

(i) Type 1 - Velocity potentials are specified: 

(|) = q^x (99) 
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Figure 24. Grid C (Fine - 40x20) 


(ii) 


Type 2 - Fluxes are specified: 


f 


a 


(100) 


The surface distribution of the pressure coefficients obtained using all 
three grids are as given in figures 25 and 26. Results of the intermediate 
and fine grids are in close agreement whereas the coarse grid results deviate 
significantly from the other two at the shock. Flux boundary conditions 
yield slightly stronger shocks and the shock positions are one element fur- 
ther downstream of what was predicted by the velocity potential boundary 
conditions . 

In order to determine the effects of replacing an unbounded domain by 
a large bounded domain, the above analyses were repeated' by taking the far- 
field boundaries of grid B also at 5R and 14R radial distances from the 
center. The grids thus obtained are shown respectively in figures 27 and 
28. As it is seen from the obtained results in figure 29, the solutions of 
lOR and 14R grids are in agreement for the velocity potential boundary con- 
dition whereas the 5R grid predicts considerable weaker shocks. The com- 
parison of the results for flux boundary conditions in figure 30 indicates 
an opposite trend and the shock in the 5R grid is stronger. 

The above results for the cylinder were obtained by employing the 

shock capturing scheme. Shock fitting applied to the radial grids shown 

in figures 22 and 24 failed. This is attributed to the misalignment of 

the element interfaces with the shock line because, for the shock fitting 

to be successful, the element interfaces must be reasonably aligned with 

the predicted shock directions. The degree of misalignment of element 

interfaces can be determined by computing the direction of the shock line 

on the surface of the cylinder using isentropic shock-jump relations. For 

instance, from the computed quantities on two sides of the shock for the 

grid shown in figure 23 with velocity potential boundary conditions, the 

shock direction at the surface of the cylinder is predicted to be 85 with 

respect to the x-axis (positive counter-clockwise) whereas the interface 

of the element at that location makes an angle of 60° with the x-axis. 

For this reason the radial grid B is modified so that the element interfaces 

in the supersonic pocket are more in line with the shock direction as shown 

in figure 31. With the new results in this case, the shock direction is 
o o 

computed to be 92 while the element interfaces make an angle of 86 which 

indicates a better alignment. The comparison of shock fitting and shock 
capturing results for this case is presented in figure 32. As it is ob- 
served, both the shock fitting and the shock capturing results of the 
modified grid are in good agreement, indicating the accuracy of the solutions. 
Moreover, slight smearing present at the upstream shock region of the cap- 
turing solution in the radial grid disappears with the modification of the 
grid. For visualization of the supersonic pocket geometry with the compu- 
tational grid, the sonic line is superimposed both for the regular and the 
modified grids as shown respectively in figures 23 and 31. Also, iso-Mach 
and iso-potential lines of this problem are shown in figures 33 and 34. 

Mach contours in figure 33 are plotted with increments of Z\M = 0.1. 
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Figure 28. Grid B with 14R outer radius (30x15) 
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Figure 32. Comparison of shock fitting and shock 
capturing solutions 
(M = 0.51, 4) specified) 







Figure 34. Iso-Potential lines for the modified grid 
(M^ = 0.51, (}) specified) 


It may be observed from figure 34 that the modified grid is better aligned 
with the constant potential lines in the supersonic pocket. 

The final results obtained for this M =0.51 farfield condition are 

oo 

summarized in figure 35. Shown in the same figure is the result reported 
by Hafez, et.al., (ref. 19) who employed a finite difference-artificial 
compressibility method. For the same problem they have used 60x30 
uniform radial grid of 5R outer radius and specified velocity potentials at 
the farfield boundaries using a farfield formula. 

As may be seen, the finite element results obtained for the various 
boundary conditions, including the result obtained with the addition of a 
source term to the uniform flow expression in equation (99), are all in 
close agreement. However, there is a considerable difference between the 
present results and the result reported in Reference 19 . 

The application of an additional source term 

<l>a = ‘Ico ^ (101) 

obtained from the incompressible potential solution for the farfield did 
not seem to affect the results of lOR and 14R grids to any noticeable degree 
but increased the strength of the shock considerably in the 5R grid. 

b) Case Ib : M =0.45. -The effect of grid refinement on the results is 

summarized in figure 36 for this flow condition. Flux boundary conditions 
and lOR radial grids were used in each case. As can be seen, both inter- 
mediate and fine grids predict essentially the same solution whereas the 
coarse grid again does not have the required accuracy due to poor resolution 
near the shock. 

Shown in figure 37 are the surface distribution of pressure coefficients 
obtained using 5R, lOR and 14R intermediate size grids. Also shown in the 
same figure for comparison is the solution reported by Bristeau, et al. , 

(ref. 6) who employed a conjugate gradient solution algorithm in conjunction 
with finite elements having a grid of 3456 triangular elements and 1813 nodal 
points. They have used flux boundary conditions on a large bounded polygonal 
domain; however, the size of the outer boundary was not indicated. As it may 
be seen, the present 5R grid solution agrees well with the solution of 
reference 6. 
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Figure 36. Effect of grid refinements 
(M =0.45, flux specified) 

OO 


74 




Test Case II 

45° Staggered Low Solidity NACA 0012 Cascade 

A cascade of NACA 0012 airfoils with a 45° stagger angle and 
h/C = 3.6 is analyzed for the following inlet flow conditions: 


Case 

Ila: 

M. =0.8 
in 

9 

3- = 1° 

in 

Case 

lib: 

M. =0.8 
in 

9 

3. = 0° 
in 

Case 

lie: 

M. =0.8 
in 

9 

3. =-l' 

in 


Kutta condition was applied at the trailing edge for the above cases. 

The chordwise distribution of the pressure coefficients obtained for 
each case using the finite element grid shown in figure 4 on page 19 are 
respectively given in figures 38-40. There are only 20 elements over 
each surface of the airfoil in this grid. The plotted results are those 
computed at the centroid of each element. With the combined shock capturing 
and shock fitting scheme, convergent unique solutions were reached at 

= 1.25, 1.5 and 4.2 respectively for Cases Ila, Ilb, and lie. For these 

cases the maximum computed local Mach numbers were respectively = 1.25, 

1.16 and 1.07 over the upper surface, and =1.11, 1.26 and 1.39 over 

the lower surface of the airfoil. 

The need for large artificial viscosity in Case lie is attributed to 
the size of the supersonic pocket, as has been implied by the von Neumann 
convergence analysis. The greater the size of the supersonic pocket, the 
higher the amount of artificial viscosity required for uniform convergence. 
This general trend is evident in all three cases analyzed herein. 

The development of sharp shock fronts with almost no smearing and the 

satisfaction of the trailing edge Kutta conditions are apparent in all 

cases. In addition, figures 38—40 show that, with the above set of 

artificial viscosity multipliers, y , more smearing is produced along the 

0 

upper surface of the airfoil regardless of the local Mach number. This may 
be attributed to the slight misalignment of the grid with the local flow 
direction. 

Iso-Mach line plots of the above three cases are presented in figures 
41-43. On these figures sonic lines are denoted as darkened contours. The 
size of the supersonic pocket and the extension of the shock beyond the 
mid-channel of the cascade in Case lie are of interest. In this case, with 
the higher values of artificial viscosity, convergent solutions were obtained 
initially with the shock placed close to the mid-chord of the airfoil. As 
the artificial viscosity was reduced, the shock moved downstream and the 
size of the supersonic pocket was increased. Each time the artificial vis- 
cosity was reduced, shock fitting was applied to the convergent shock 
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Figure 39. 


Variation of pressure coefficients 

the airfoil for Case Ilb (M. = 0 

in 


3. = 0°, M 

in upper max 


1.16, M 


over 

8 , 


lower max 



Figure 40. Variation of pressure coefficients 
over the airfoil for Case lie 

(M. = 0.8, 3. = -1°, M 

in upper max 

M. 

lower max 


= 1.39) 
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capturing solution. This was repeated until the shock moved downstream 
far enough to satisfy the shock jump conditions; in which case both the 
shock fitting and the shock capturing yielded the same solution. 

Finally, the chordwise pressure distributions obtained for an 
additional test case are presented for upper and lower surface of the 
airfoil respectively in figures A4 and 45. Here the inlet flow conditions 
of Case lib were used, but instead of determining the exit flow angle from 
the trailing edge Kutta condition, a value of 3^^ 0 was specified along 

the exit stations. Both the grid shown in figure 4 (20 elements over the 
airfoil) and a fine grid shown in figure 46 (30 elements over the airfoil) 
were employed for this analysis. The differences in the pressure distri- 
butions obtained from these grids are too small to show in figures 44 and 45 
Also shown in the same figures, are the results of a finite difference- arti 
ficial compressibility scheme employed by Farrell and Adamczyk (ref. 16). 

The results obtained by both methods are within modeling accuracy. 


Test Case III 

45° Staggered High Solidity NACA 0012 Cascade 

A high solidity cascade of NACA 0012 airfoils with a stagger angle of 
45° and h/C =1, is considered in this case. The computational grid em- 
ployed in the analysis is given in figure 47. 

The flow through this cascade is analyzed for the following inlet 
Mach numbers and angles of attack: 


Case 

Ilia: 

M. =0.80, 
in 

3. = 0°, 

in 

3 = 4.01° 

ex 

(Computed) 

Case 

mb: 

M. = 0.85, 
in 

= 0°. 

3 = 4.70° 

ex 

(Computed) 

Case 

IIIc: 

M. =0.85, 
in 


3 = 4.60° 

ex 

(Computed) 

Case 

Illd: 

M. = 0.87, 
in 

^in = 

3 = 4.87° 

ex 

(Computed) 


The obtained numerical results for the pressure coefficients and iso-Mach 
lines are presented in figures 48 and 49 respectively. Mach contours in 
figure 49 are plotted with increments of AM = 0.05 and the sonic lines on 
these figures are denoted by darkened lines. 
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Figure 45. Variation of pressure coefficients over the 
lower surface on the airfoil for a specified 

exit angle (M. = 0.8, 3- =6 = 0^) 

in in ex ' 


Figure 46. Refined grid employed for Case II with a specified exit angle 
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Figure 47. Finite element grid for 45° staggered, high solidity cascade 
of NACA 0012 


a) Case Ilia (M, = 0.80, 3. - 0 ) 
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b) Case Illb (M. = 0.85, 3 
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Figure 48. Pressure coefficients over the airfoil surfaces for Case III 



c) Case IIIc (M. = 0.85, 3. = 1 ) 

in in 


d) Case Illd (M. = 0.87, B 

in 


Figure 48. Pressure coefficient over the airfoil surfaces for Case III (Continued) 



a) Case Ilia (M. = 0.80, 3. = 0 ) b) Case Illb (K. = 0.85, 3 

in in in 

-Choked flow 


Figure 49. Iso-Mach lines for Case III ( AM = 0.5 ) 



a) Case IIIc (M. = 0.85, S = 1°) 

in in 

-Choked flow 


Figure 49. Iso-Mach lines for Case III (Continued) 




For Case Ilia, the solution consists of embedded shocks, one on each 
surface of the airfoil, the lower shock being much stronger. As the inlet 
Mach number is increased to 0.85 in Case Illb, the strength and the size 
of both shocks grow considerably. Meanwhile, a new supersonic pocket de- 
velops on the upper surface. At this point, a subsonic to supersonic sonic 
line crosses the entire channel and the flow is choked. It is interesting 
to note the increased turning effect of the Kutta condition at the trailing 
edge as seen in figure 49b . 

When the inlet Mach number is further increased to 0.87, a new set of 
stable solutions are obtained. As can be observed from figures 49b, 49d, 
and 50, the position of the sonic line crossing the channel does not change. 
In fact, the flow field past the sonic line in the channel, remains virtu- 
ally the same except in the vicinity of the shocks. Since, the outlet mass 
flow is increased for M^^ = 0.87, the flow turns to produce an oblique 

shock to satisfy the exit flow conditions. The change in circulation, in 
turn, affects the flow configuration around the trailing edge. The sen- 
sitivity of flow conditions around the trailing edge will further be 
discussed . 

In the above test cases, no additional difficulties were encountered 
due to the developed complex shock patterns. The obtained results illustrate 
the generality of the computational method. Multiple shocks, supersonic 
pockets with irregular geometries, interaction of the trailing edge con- 
ditions with the shocks, and the effect of large stagger angles can all be 
accounted for. In all cases, the accuracy of the shock capturing solutions 
were tested by comparing each with the shock fitting results (ref. 15) and 
found to be satisfactory. Furthermore, it should be noted that the Mach 
contour patterns obtained for the choked flow cases are in close agreement 
with those reported by Kozel, et.al. (ref. 27). Their results include both 
finite difference solutions of small disturbance equations and experimental 
measurements for a cascade of 8%, symmetric, circular-arc airfoils with 

a = 45° and h/C = 1. 
s 


Test Case IV - 8% Circular-Arc Cascade 

Using the grid shown in figure 51, the 45° staggered cascade of 8% 

circular-arc airfoil of Kozel, et.al. (ref. 27) is also analyzed here for 

the inlet flow conditions M. = 0.876 and B. = 0.5°. Shown in figure 52 

in in 

are the experimental and numerical results reported in reference 27, as 
well as the finite element results of the same problem. The agreement of 
the numerically determined iso-Mach lines with that of the experimentally 
obtained iso-density lines is excellent. The calculated pressure distri- 
bution over the airfoil surface is shown in figure 53. The exit flow angle 
for this case was computed to be 2.88°. 
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Figure 51. 


Finite element grid for 45 


o 


cascade of 8% circular-arc airfoil 
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Figure 52. Iso-Mach and iso-density lines for 8% 
circular-arc airfoil 

a) experimental (ref. 27), 

b) finite difference solution of small 
perturbation equation (ref. 27), 

c) finite element solution ( AM = 0.5 ) 





Test Case V - Supercritical Cascade 

Shown in figure 54 is the computational grid used for a cascade of 
shock-free supercritical blades designed by Sanz (ref. 30) from the hodo- 
graph solution. The cascade was analyzed using the following shock-free 
design conditions: 

M. =0.711 
xn 

3 . =30.81° 

xn 


a =9.33° 
s 

h/C = 1.034 . 

The same cascade was also analyzed by Steger, et.al., (ref. 32) using an 
implicit finite difference code for full-Euler equations. 

Mach number distributions obtained from the three methods are 
presented in figure 55 for comparison, showing the present solution to be 
in excellent agreement with the hodograph solution. The Mach number con- 
tours of the flow, determined by the present finite element method, are 
displayed in figure 56 using AM = 0.05 increments. 


No particular difficulty was encountered in obtaining steady state 
solutions to this problem. A small "wiggle” appears in the finite element 
solution of the pressure distribution at the upper surface of the blade 
which does not exist in the hodograph solution. This is attributed to the 
sensitivity of the shock-free solution to the accurate definition of the 
geometry at this point of high curvature. A refinement of the grid at this 
region seems necessary for better accuracy. 


Efficiency of the Developed Computational Procedure 


The solution of the transonic flow problem in a cascade by using the 
developed procedure involves the following basic steps: 

. generate a finite element grid, 

. start with a high artificial viscosity parameter, 
iterate until a steady state solution is obtained, 

. reduce the artificial viscosity and repeat the procedure. 
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Figure 54. Computational grid used for a cascade of shock-free 
blade of Sanz 
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. perforin shock fitting when "wiggles" occur at the 
shock region, 

. repeat the procedure with a more refined grid to check 
the accuracy of the obtained solution if necessary. 

The computational efficiency for each of the above steps are based on 
performing the following tasks: 

. the generation of a finite element grid requires the 
definition of nodal coordinates for each node and nodal 
connectivities for each element. In the applications 
summarized in this report, the grid was prepared by coding 
a simple special purpose mesh generation program for each 
case. In practice, this can be achieved interactively by 
using general purpose finite element mesh generation 
programs and computer graphics terminals, e.g., (ref. 34), 

. iterative solution of the equations by using the constant 
coefficient scheme described in page 13 requires decom- 
position of a symetric matrix only once at the beginning 
of the iterations. This matrix is banded for which a 
wave-front solver is employed. After the first time step, 
back-substitution procedure is performed on this system 
of equations, 

. the computer cost is proportional to the number of nodal 
points and elements on the computational grid, and the 
number of iterations required for final convergence. These 
depend on the complexity of the problem and the users 
knowledge of the problem. For example when one solves a 
problem for 0 angle of attack for a cascade, the solution 
of 1 angle of attack case can be obtained by choosing 
better initial conditions and a more efficient variation 
of artificial viscosity and relaxation factor during the 
time steps. 

The computer program named TRANS4 which has evolved during this re- 
search has been coded and tested for both CDC 6600 and IBM 4341 machines. 
The central and auxiliary memory requirements of the program, in number 
of words, for CDC 6600 can be determined by using the expressions: 

CM = Max {8*N0DES, 6*N0DES + WF*(WF-l)/2 + 10 } + 5*NUMEL + 35,500 

MS = 53-(NUMEL + NODES) 

where 

NODES = total number of nodal points, 

NUMEL = total number of elements. 
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WF = maximum wave-front of equations, 

CM = central memory requirement, 

MS = information to be stored on auxiliary memory devices. 

In order to give a detailed description of the computer cost, computer 
requirements of four sample cases are presented in Table 1. The total 
CPU times are quoted based on the assumptioij that uniform flow conditions 
are used as initial solutions in each case and that no additional flow 
information is available. 
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CONCLUSIONS 


Application of the finite element method to the solution of transonic 
flow problem in cascades is presented and the importance of the computa- 
tional grid for obtaining accurate and efficient solutions is illustrated. 
The following basic conclusions, which were reached through this study, 
summarize the importance considerations in the solution of transonic flow 
problems in cascades: 

. The most critical regions in the modeling of transonic flows 
are the supersonic pocket and the shock. The accuracy and 
convergence on the solution in this region primarily depend 
on the artificial viscosity distribution. The artificial 
viscosity distribution must be chosen in terms of the local 
Mach number, the local grid size, and the number of elements 
in the supersonic pocket. 

. The shock is a critical region since errors in the supersonic 
pocket are magnified as they are convected in the direction 
of the flow till they reach the shock. The solution is most 
sensitive to the artificial viscosity distribution at this 
region. A combined usage of shock fitting and shock capturing 
techniques is very useful for the solution of the problem. 
Finite element method is very efficient in applying these 
techniques . 

. The leading edge of the airfoil requires a refined grid for 
better accuracy. However, if the flow is subsonic, the 
smearing of the solution in this region with a coarse grid 
remains local and does not generally affect the overall flow 
configuration, since the equations are elliptic in this 
region. 

. Application of the Kutta condition illustrates the sensitivity 
of the flow around the trailing edge region and the importance 
of proper exit conditions. The need for accurate modeling on 
the trailing edge is apparent in the sample cases in which 
interaction with the supersonic pocket develops. 

. Development of a finite element grid for the problem requires 
consideration of several different flow structures. If one 
uniformly refines the grid for accuracy, the efficiency suffers 
considerably, since the convergence requirements at the super- 
sonic pocket become more stringent with the increasing number 
of elements. Since the pressure gradients are generally low, 
one should stay with a reasonable number of elements in this 
region to obtain accurate yet efficient results. Leading and 
trailing edges are local regions of high pressure gradients and 
require local refinements in the grid. The modeling of the 
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shock requires a grid fitting rather than a grid refinement. 

This suggests a modification on the grid rather than a refine- 
ment . 

The obtained numerical results demonstrate that the finite element 
method can be employed quite effectively for the analysis of rather complex 
transonic flow problems in cascades. Simple grids that will provide accurate 
results as well as fast convergence can be generated quite efficiently. The 
above discussion also describes the basic advantage of a finite element 
method for the solution of transonic flow problems. In addition to the compu- 
tat"ional grids which fit the boundaries of the overall flow domain, one 
should also be able to generate grids which are capably of accurately 
describing the irregi.ilar geometries of different flow regions. The grid 
generation should be interactive to allow modifications at critical regions, 
such as around the shock or at the leading and trailing edges. The current 
finite element mesh generation techniques provide such capabilities resulting 
with considerable computational efficiencies in the solution of transonic 
flow problems in cascades, e.g., (ref. 34). 

Finally, for the designer of a new cascade, the following analysis 
procedure can be recommended: 

. Start with a coarse grid and calculate a basic solution 
to the problem. 

. Check the accuracy of the solution by refining the grid at 
critical regions and comparing the results. 

Once an efficient grid is defined for a particular set of problems, 
the designer can approach the problem with more confidence. In terms of 
practical efficiency in the solution of transonic flow problems, i.e., in 
terms of the number of man-hours and amount of computer time required to 
obtain an accurate solution to the problem, the above approach is practical. 
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